#pragma once
#include <cuda_runtime.h>
#include <glm/glm.hpp>

namespace PhysLeo {

/**
 * peridynamic anisotropic nonlinear material.
 * \f$ g=h=\sum^4_{i=1}{p_iA_i} , G=|R\sum v|\f$.
 * \f$A_i\f$ is A class nonlinear basis function, \f$p_i\f$ is the nonlinear parameters which controls the combination of different nonlinear basis.
 * R is a rotation matrix, \f$\sum\f$ is a diagnal matrix, v is original bond direction.
 */
template<typename T>
class PdmAnisoNonlinear {
public:
    /**
    * g function in Pdm anisotropic nonlinear material is \f$ \sum^4_{i=1}{p_iA_i} \f$
    * @param[in] x  a T value describe the deformation severity
    * @param[in] ptr_material  a pointer to vertex specific material data
    * @return T value
    */
    __host__ __device__ static T gFunction(T x, T* ptr_material)
    {
        T nonlinear_params[4] = { ptr_material[2],ptr_material[3] ,ptr_material[4] ,ptr_material[5] };

        auto x2 = x * x;
        auto x3 = x2 * x;
        auto x4 = x2 * x2;
        auto x5 = x4 * x;
        auto sum = nonlinear_params[0] * (0.5f *x2 - log(x) - 0.5f);
        sum += 0.5f*nonlinear_params[1] * (1.0f / 3.0f*x3 + 1.0f / x - 4.0f / 3.0f);
        sum += 1.0f / 3.0f*nonlinear_params[2] * (1.0f / 4.0f*x4 + 1.0f / 2.0f / x2 - 3.0f / 4.0f);
        sum += 1.0f / 4.0f*nonlinear_params[3] * (1.0f / 5.0f*x5 + 1.0f / 3.0f / x3 - 8.0f / 15.0f);
        return sum;
    }

    /**
    * h function in Pdm anisotropic nonlinear material is \f$ \sum^4_{i=1}{p_iA_i} \f$
    * @param[in] x  a T value describe the deformation severity
    * @param[in] ptr_material  a pointer to vertex specific material data
    * @return T value
    */
    __host__ __device__ static T hFunction(T x, T* ptr_material)
    {
        return gFunction(x, ptr_material);
    }

    /**
    * first derivative of g function in Pdm anisotropic nonlinear material is \f$ \sum^4_{i=1}{p_iA_i^{'}} \f$
    * @param[in] x  a T value describe the deformation severity
    * @param[in] ptr_material  a pointer to vertex specific material data
    * @return T value
    */
    __host__ __device__ static T dgFunction(T x, T* ptr_material)
    {
        T nonlinear_params[4] = { ptr_material[2],ptr_material[3] ,ptr_material[4] ,ptr_material[5] };

        auto x2 = x * x;
        auto x3 = x2 * x;
        auto x4 = x2 * x2;

        auto sum = nonlinear_params[0] * (x - 1.0f / x);
        sum += 0.5f*nonlinear_params[1] * (x2 - 1.0f / x2);
        sum += 1.0f / 3.0f*nonlinear_params[2] * (x3 - 1.0f / x3);
        sum += 0.25f*nonlinear_params[3] * (x4 - 1.0f / x4);
        return sum;
    }

    /**
    * first derivative of h function in Pdm anisotropic nonlinear material is \f$ \sum^4_{i=1}{p_iA_i^{'}} \f$
    * @param[in] x  a T value describe the deformation severity
    * @param[in] ptr_material  a pointer to vertex specific material data
    * @return T value
    */
    __host__ __device__ static T dhFunction(T x, T* ptr_material)
    {
        return dgFunction(x, ptr_material);
    }

    /**
    * second derivative of g function in Pdm anisotropic nonlinear material is \f$ \sum^4_{i=1}{p_iA_i^{''}} \f$
    * @param[in] x  a T value describe the deformation severity
    * @param[in] ptr_material  a pointer to vertex specific material data
    * @return T value
    */
    __host__ __device__ static T ddgFunction(T x, T* ptr_material)
    {
        T nonlinear_params[4] = { ptr_material[2],ptr_material[3] ,ptr_material[4] ,ptr_material[5] };

        auto x2 = x * x;
        auto x3 = x2 * x;
        auto x4 = x2 * x2;
        auto x5 = x4 * x;
        auto sum = nonlinear_params[0] * (1.0f + 1.0f / x2);
        sum += nonlinear_params[1] * (x + 1.0f / x3);
        sum += nonlinear_params[2] * (x2 + 1.0f / x4);
        sum += nonlinear_params[3] * (x3 + 1.0f / x5);

        return sum;
    }

    /**
    * second derivative of h function in Pdm anisotropic nonlinear materialn is \f$ \sum^4_{i=1}{p_iA_i^{''}} \f$
    * @param[in] x  a T value describe the deformation severity
    * @param[in] ptr_material  a pointer to vertex specific material data
    * @return T value
    */
    __host__ __device__ static T ddhFunction(T x, T* ptr_material)
    {
        return ddgFunction(x, ptr_material);
    }

    /**
    * anisotropic weight function in Pdm anisotropic nonlinear material, |Gv|.
    * G is the 3x3 anisotropic matrix, it can be construct by \f$ R\sum \f$, R is rotation matrix represents three main direction, \f$ \sum \f$ is diagonal matrix represents the specific intensity in different direction.
    * @param[in] direction  a vec3<T> value describe the material direction
    * @param[in] ptr_material  a pointer to the specific vertex material data
    * @return T value
    */
    __host__ __device__ static T anisoWeight(glm::tvec3<T> direction, T* ptr_material)
    {
        auto aniso = glm::tmat3x3<T>(
            ptr_material[6], ptr_material[7], ptr_material[8],
            ptr_material[9], ptr_material[10], ptr_material[11],
            ptr_material[12], ptr_material[13], ptr_material[14]);

        return glm::length(aniso*direction);
    }
};

}